# Import Pydicom before using Pydicom functions
# DICOM (Digital Imaging and Communications in Medicine) is the standard protocol for
# the management and transmission of medical images and related data, and is used by
# many healthcare facilities
# Pydicom is a pure Python package for working with DICOM files
import pydicom
# The `dicom2nifti` library converts DICOM files to the NIfTI file and supports most anatomical
# CT and MRI data, especially MRI, which supports most 4D data such as DTI data and fMRI data
import dicom2nifti
# NiBabel is a pure Python package that supports a growing number of neuroimaging file formats,
# including the NIfTI1 format and the NIfTI2 format, which is an extension of the former
# NiBabel provides both format-independent access to high-level neuroimaging and format-specific APIs with
# different levels of access to all available information in a particular file format
# NiBabel's API provides full or selective access to file header information (metadata), while image data
# is provided via NumPy arrays
import nibabel as nib
# The `pathlib` module is similar to the `os.path` module, but `pathlib` provides a more
# advanced and convenient interface than `os.path`
# It is possible to use `pathlib` to represent file paths as specialized `Path` objects
# instead of plain strings
from pathlib import Path
# SimpleITK is a simplified programming interface to the algorithms and data structures
# of the Insight Toolkit (ITK) for segmentation, registration and advanced image analysis
# SimpleITK supports bindings for multiple programming languages including C++, Python, R,
# Java, C#, Lua, Ruby and TCL
# Combining SimpleITK’s Python bindings with the Jupyter notebook web application creates
# an environment which facilitates collaborative development of biomedical image analysis
# workflows
# In this project, SimpleITK provides the ability to automatically detect and read all DICOM
# files so that the user does not have to manage file reading or slice sorting
import SimpleITK as sitk
import numpy as np
import torch
from torchvision.utils import make_grid
import matplotlib as mpl
import matplotlib.gridspec as gridspec
import matplotlib.pyplot as plt
from datetime import datetime
from functools import wraps
import itertools
import math
import random
import reprlib
from scipy import ndimage
import sys
%matplotlib inline
XINHUI = "#7a7374"
XUEBAI = "#fffef9"
YINBAI = "#f1f0ed"
YINHUI = "#918072"
figure_size = (16, 9)
custom_params = {
"axes.axisbelow": True,
"axes.edgecolor": YINBAI,
"axes.facecolor": XUEBAI,
"axes.grid": True,
"axes.labelcolor": XINHUI,
"axes.spines.right": False,
"axes.spines.top": False,
"axes.titlecolor": XINHUI,
"figure.edgecolor": YINBAI,
"figure.facecolor": XUEBAI,
"grid.alpha": 0.8,
"grid.color": YINBAI,
"grid.linestyle": "--",
"grid.linewidth": 1.2,
"legend.edgecolor": YINHUI,
"patch.edgecolor": XUEBAI,
"patch.force_edgecolor": True,
"text.color": XINHUI,
"xtick.color": YINHUI,
"ytick.color": YINHUI,
}
mpl.rcParams.update(custom_params)
reprlib_rules = reprlib.Repr()
reprlib_rules.maxother = 250
sys.path.append("../")
from Modules import *
# The dataset used in this practice project is a very small subset of CT images extracted from
# the Cancer Imaging Archive (TCIA), which contains intermediate slices of all CT images
# in which valid age, mode, and contrast markers can be found
# This dataset is provided by the Kaggel Datasets called CT Medical Imaging
# (https://www.kaggle.com/datasets/kmader/siim-medical-images),
# license type is Attribution 3.0 Unported (CC BY 3.0)
# https://creativecommons.org/licenses/by/3.0/legalcode
dir_path = "../Datasets/Kaggle - CT Medical Images/dicom_dir/"
sample_dcm = "ID_0000_AGE_0060_CONTRAST_1_CT.dcm"
# The main function of Pydicom to read and parse DICOM files is `read_file`
dicom_file = pydicom.read_file(dir_path + sample_dcm)
tabulation = Form_Generator()
tabulation.heading_printer("Initial understanding of DICOM file")
statements = [
"""
dir_path = "../Datasets/Kaggle - CT Medical Images/dicom_dir/"
sample_dcm = "ID_0000_AGE_0060_CONTRAST_1_CT.dcm"
dicom_file = pydicom.read_file(dir_path + sample_dcm)
"""
]
tabulation.statement_generator(statements)
variables = ["dicom_file"]
values = [str(reprlib_rules.repr(dicom_file))]
tabulation.variable_generator(variables, values, 1)
expressions = [
"dicom_file[0x0028, 0x0010]",
"dicom_file[0x0028, 0x0011]",
"dicom_file[0x0018, 0x0015]",
"dicom_file.Rows",
"dicom_file.Columns",
"dicom_file.BodyPartExamined",
"dicom_file.keys()",
"dicom_file.values()",
"dicom_file.dir()",
"dicom_file.dir('Image')",
]
results = [
str(dicom_file[0x0028, 0x0010]),
str(dicom_file[0x0028, 0x0011]),
str(dicom_file[0x0018, 0x0015]),
str(dicom_file.Rows),
str(dicom_file.Columns),
str(dicom_file.BodyPartExamined),
str(reprlib_rules.repr(dicom_file.keys())),
str(reprlib_rules.repr(dicom_file.values())),
str(reprlib_rules.repr(dicom_file.dir())),
str(reprlib_rules.repr(dicom_file.dir("Image"))),
]
tabulation.expression_generator(expressions, results, 1)
Initial understanding of DICOM file +-------------------------------------------------------+ | Statement | +-------------------------------------------------------+ | dir_path = "../Datasets/Kaggle - CT Medical | | Images/dicom_dir/" | | sample_dcm = "ID_0000_AGE_0060_CONTRAST_1_CT.dcm" | | dicom_file = pydicom.read_file(dir_path + sample_dcm) | +-------------------------------------------------------+ +------------+------------------------------------------------+ | Variable | Value | +------------+------------------------------------------------+ | dicom_file | Dataset.file_meta | | | ------------------------------- | | | (0002, 0000) File Meta Information Group | | | Length UL: 194 | | | (0002, 0001) Fil...Group Length | | | UL: 524296 | | | (7fe0, 0010) Pixel Data | | | OW: Array of 524288 elements | +------------+------------------------------------------------+ +-----------------------------+-------------------------------+ | Expression | Result | +-----------------------------+-------------------------------+ | dicom_file[0x0028, 0x0010] | (0028, 0010) Rows | | | US: 512 | | dicom_file[0x0028, 0x0011] | (0028, 0011) Columns | | | US: 512 | | dicom_file[0x0018, 0x0015] | (0018, 0015) Body Part | | | Examined | | | CS: 'CHEST' | | dicom_file.Rows | 512 | | dicom_file.Columns | 512 | | dicom_file.BodyPartExamined | CHEST | | dicom_file.keys() | dict_keys([(0008, 0000), | | | (0008, 0005), (0008, 0008), | | | (0008, 0016), (0008, 0018), | | | (0008, 0020), (0008, 0021), | | | (0008, 0022), ...095, 0010), | | | (0097, 0000), (0097, 0010), | | | (0099, 0000), (0099, 0010), | | | (7003, 0000), (7003, 0010), | | | (7fe0, 0000), (7fe0, 0010)]) | | dicom_file.values() | dict_values([(0008, 0000) | | | Group Length | | | UL: 430, (0008, 0005) | | | Specific Character Set | | | CS:...up Length | | | UL: 524296, | | | (7fe0, 0010) Pixel Data | | | OW: | | | Array of 524288 elements]) | | dicom_file.dir() | ['AccessionNumber', | | | 'AcquisitionDate', | | | 'AcquisitionNumber', | | | 'AcquisitionTime', | | | 'BitsAllocated', | | | 'BitsStored', ...] | | dicom_file.dir('Image') | ['ImageDimensions', | | | 'ImageFormat', | | | 'ImageGeometryType', | | | 'ImageLocation', | | | 'ImageOrientation', | | | 'ImageOrientationPatient', | | | ...] | +-----------------------------+-------------------------------+
# DICOM data has an attribute called `pixel_array` that provides more useful pixel data
# for uncompressed images that can be passed to the graphics library for viewing
# To use this attribute, the system must have the NumPy numeric package installed,
# since `pixel_array` returns a NumPy array
ct = dicom_file.pixel_array
tabulation = Form_Generator()
tabulation.heading_printer("Getting pixel data from DICOM file")
statements = ["ct = dicom_file.pixel_array"]
tabulation.statement_generator(statements)
variables = ["ct"]
values = [str(reprlib_rules.repr(ct))]
tabulation.variable_generator(variables, values)
expressions = ["ct.shape"]
results = [str(ct.shape)]
tabulation.expression_generator(expressions, results)
Getting pixel data from DICOM file +-----------------------------+ | Statement | +-----------------------------+ | ct = dicom_file.pixel_array | +-----------------------------+ +----------+------------------------------------------------+ | Variable | Value | +----------+------------------------------------------------+ | ct | array([[0, 0, 0, ..., 0, 0, 0], | | | [0, 0, 0, ..., 0, 0, 0], | | | [0, 0, 0, ..., 0, 0, 0], | | | ..., | | | [0, 0, 0, ..., 0, 0, 0], | | | [0, 0, 0, ..., 0, 0, 0], | | | [0, 0, 0, ..., 0, 0, 0]], dtype=uint16) | +----------+------------------------------------------------+ +------------+------------+ | Expression | Result | +------------+------------+ | ct.shape | (512, 512) | +------------+------------+
def image_display(image, ax, title, cmap):
ax.imshow(image, cmap)
ax.grid(False)
ax.set_title(title, loc="center", pad=10)
x_ticks = list(range(0, image.shape[1], 50))
y_ticks = list(range(0, image.shape[0], 50))
ax.set(xticks=x_ticks, xticklabels=x_ticks,
yticks=y_ticks, yticklabels=y_ticks)
ax.set_xlim(left=0)
ax.set_ylim(top=0)
ax.minorticks_on()
return ax
# Path classes are divided between pure paths, which provide purely computational
# operations without I/O, and concrete paths, which inherit from pure paths but also
# provide I/O operations
path_object = Path(dir_path)
# `PurePath.name` returns a string representing the final path component, excluding
# the drive and root directory (if any)
# When a path points to a directory, `Path.iterdir` generates a path object of the directory's
# contents
# `Path.is_file` returns True if the path points to a normal file or to a symbolic link
# to a normal file, False if it points to another type of file
random_dicom_path = random.choice(
[
file.name
for file in path_object.iterdir()
if file.is_file() & (pydicom.read_file(file).BodyPartExamined != "CHEST")
]
)
random_dicom_file = pydicom.read_file(dir_path + random_dicom_path)
plt.rcParams["figure.figsize"] = (
figure_size[0] / 3 * 2, figure_size[1] / 4 * 5)
fig, axs = plt.subplots(nrows=2, ncols=2)
image_display(
ct,
axs[0, 0],
f"{dicom_file.BodyPartExamined.capitalize()} CT medical image displayed using "
"colormap 'gray'",
cmap="gray",
)
image_display(
ct,
axs[0, 1],
f"{dicom_file.BodyPartExamined.capitalize()} CT medical image displayed using "
"colormap 'bone'",
cmap="bone",
)
image_display(
random_dicom_file.pixel_array,
axs[1, 0],
f"{random_dicom_file.BodyPartExamined.capitalize()} CT medical image displayed using "
"colormap 'gray'",
cmap="gray",
)
image_display(
random_dicom_file.pixel_array,
axs[1, 1],
f"{random_dicom_file.BodyPartExamined.capitalize()} CT medical image displayed using "
"colormap 'bone'",
cmap="bone",
)
fig.suptitle(
"Visual Comparison of CT Medical Image Display Using Different Colormaps",
fontsize="x-large",
x=0.5,
y=0,
)
plt.tight_layout()
plt.show()
def annotation_color(cmap):
cmap_color = mpl.colormaps[cmap]
return cmap_color(0.85)
def first_max(list):
max_element = list[0]
for element in list[1:]:
if len(element) > len(max_element):
max_element = element
return max_element
def text_aligner(text):
line_list = text.split("\n")
max_length = len(first_max(line_list))
line_list = [
(":" + " " * (max_length - len(line))).join(line.split(":"))
if len(line) < max_length
else line
for line in line_list
]
text = "\n".join(line_list)
return text
def plane_information_annotator(func):
@wraps(func)
def wrapper(ax, cmap, dicom_file, fontsize_adjustment):
plane = func(ax, cmap, dicom_file, fontsize_adjustment)
text = f"Image: {dicom_file.Columns} x {dicom_file.Rows}"
text += f"\nImage Plane: {plane}"
text += f"\nPatient Position: {dicom_file.PatientPosition}"
text += f"\nSlice Thickness: {dicom_file.SliceThickness:.1f} mm"
text += f"\nSlice Location: {dicom_file.SliceLocation:.1f} mm"
ax.text(
0.575,
0.825,
text_aligner(text),
transform=ax.transAxes,
fontsize=9 + fontsize_adjustment,
fontfamily="monospace",
c=annotation_color(cmap),
)
return wrapper
def orientation_annotator(func):
@wraps(func)
def wrapper(ax, cmap, dicom_file, fontsize_adjustment):
plane, top, right, bottom, left = func(dicom_file)
ax.text(
0.4875,
0.95,
top,
transform=ax.transAxes,
fontweight="heavy",
fontsize=11 + fontsize_adjustment,
fontfamily="monospace",
c=annotation_color(cmap),
)
ax.text(
0.4875,
0.025,
bottom,
transform=ax.transAxes,
fontweight="heavy",
fontsize=11 + fontsize_adjustment,
fontfamily="monospace",
c=annotation_color(cmap),
)
ax.text(
0.95,
0.4875,
right,
transform=ax.transAxes,
fontweight="heavy",
fontsize=11 + fontsize_adjustment,
fontfamily="monospace",
c=annotation_color(cmap),
)
ax.text(
0.025,
0.4875,
left,
transform=ax.transAxes,
fontweight="heavy",
fontsize=11 + fontsize_adjustment,
fontfamily="monospace",
c=annotation_color(cmap),
)
return plane
return wrapper
def image_orientation(func):
@wraps(func)
def wrapper(*args, **kwargs):
Left, Right, Anterior, Posterior, Head, Feet = "L", "R", "A", "P", "H", "F"
axial_plane = [Anterior, Right, Posterior, Left]
coronal_plane = [Head, Left, Feet, Right]
sagittal_plane = [Head, Anterior, Feet, Posterior]
plane = func(*args, **kwargs)
# Patient Position (0018,5100) specifies the position of the patient relative to
# the space of the imaging device and is used for annotation purposes only,
# but does not provide a precise mathematical relationship between the patient and
# the imaging device
patient_position = dicom_file.PatientPosition
if plane == "Axial":
term_list = ["S", "DR", "P", "DL"]
top = 4 - term_list.index(patient_position[2:])
if patient_position.startswith("HF"):
right = top + 1
left = right + 2
else:
left = top + 1
right = left + 2
bottom = top + 2
top, right, bottom, left = [
axial_plane[index] if index < 4 else axial_plane[index - 4]
for index in [top, right, bottom, left]
]
elif plane == "Coronal":
term_list = ["HF", "RF", "FF", "LF"]
top = 4 - term_list.index(patient_position[:2])
if patient_position.endswith("S"):
right = top + 1
left = right + 2
else:
left = top + 1
right = left + 2
bottom = top + 2
top, right, bottom, left = [
coronal_plane[index] if index < 4 else coronal_plane[index - 4]
for index in [top, right, bottom, left]
]
elif plane == "Sagittal":
term_list = ["HF", "PF", "FF", "AF"]
top = 4 - term_list.index(patient_position[:2])
if patient_position.endswith("DL"):
right = top + 1
left = right + 2
else:
left = top + 1
right = left + 2
bottom = top + 2
top, right, bottom, left = [
sagittal_plane[index] if index < 4 else sagittal_plane[index - 4]
for index in [top, right, bottom, left]
]
else:
top, right, bottom, left = np.repeat("", 4)
return plane, top, right, bottom, left
return wrapper
@plane_information_annotator
@orientation_annotator
@image_orientation
def get_plane_information(dicom_file):
# The DICOM standard defines a patient-oriented reference coordinate system (RCS) that
# enables the user to measure the position and orientation of the image relative to
# the patient
# Image Orientation (Patient) (0020,0037) specifies the cosine of the orientation of
# the first row and column relative to the patient, which should be provided as a pair,
# where the row values for the x, y, and z axes are followed by the column values for
# the x, y, and z axes
# The orientation of the axes is determined solely by the orientation of the patient, i.e.,
# the RCS and Image Orientation (Patient) (0020,0037) values specify the orientation of
# the image frame rows and columns
IOP = dicom_file.ImageOrientationPatient
row_xyz = IOP[:3]
column_xyz = IOP[3:]
# `numpy.cross` returns the cross product of two vector arrays
# The geometric interpretation of the cross product is that two vectors determine a plane,
# and the cross product points in a direction different from both vectors
# The Patient-Based Coordinate System is a right-handed system, i.e., the vector
# cross product of a unit vector along the positive x-axis and a unit vector along
# the positive y-axis is equal to a unit vector along the positive z-axis
plane = [abs(round(x)) for x in np.cross(row_xyz, column_xyz)]
# The Image Type (0008,0008) identifies important image identification features
# and it returns the Pixel Data Characteristics, Patient Examination Characteristics,
# Modality Specific Characteristics, and Implementation Specific Identifiers
# While the third value (Modality Specific Characteristics) should identify any
# specialization specific to the Image Information Object Definition (IOD), it
# and the values that follow it are not mandatory, unlike the first two values,
# and therefore some DICOM data does not have this value
if plane[0] == 1 and plane[1] == 0 and plane[2] == 0:
return "Sagittal"
elif plane[0] == 0 and plane[1] == 1 and plane[2] == 0:
return "Coronal"
elif plane[0] == 0 and plane[1] == 0 and plane[2] == 1:
return "Axial"
else:
"Unknown"
def get_patient_information(ax, cmap, dicom_file, fontsize_adjustment, patien_age):
text = f"Patient ID: {dicom_file.PatientID}"
text += f"\nPatient's Sex: {dicom_file.PatientSex}"
try:
patien_age = dicom_file.PatientAge.strip(str(0))[:-1]
except AttributeError:
patien_age = patien_age
text += f"\nPatient's Age: {patien_age} y"
text += f"\nModality: {dicom_file.Modality}"
text += f"\nBody Part Examined: {dicom_file.BodyPartExamined}"
ax.text(
0.04,
0.825,
text_aligner(text),
transform=ax.transAxes,
fontsize=9 + fontsize_adjustment,
fontfamily="monospace",
c=annotation_color(cmap),
)
def get_date_information(ax, cmap, dicom_file, fontsize_adjustment):
date = datetime.strptime(dicom_file.StudyDate, "%Y%m%d").strftime("%Y.%m.%d")
text = f"Study Date: {date}"
ax.text(
0.04,
0.05,
text,
transform=ax.transAxes,
fontsize=7 + fontsize_adjustment,
fontfamily="monospace",
c=annotation_color(cmap),
)
def get_manufacturer_information(ax, cmap, dicom_file, fontsize_adjustment):
text = f"Manufacturer: {dicom_file.Manufacturer}"
ax.text(
0.525,
0.05,
text.rjust(36),
transform=ax.transAxes,
fontsize=7 + fontsize_adjustment,
fontfamily="monospace",
c=annotation_color(cmap),
)
def DICOM_2D_image(ax, cmap, dicom_file, title, fontsize_adjustment=0, patien_age="__"):
ax.imshow(dicom_file.pixel_array, cmap)
ax.grid(False)
ax.set_title(title, loc="center", pad=10, fontsize=12 + fontsize_adjustment)
ax.set(xticks=[], yticks=[], frame_on=False)
get_patient_information(ax, cmap, dicom_file, fontsize_adjustment, patien_age)
get_plane_information(ax, cmap, dicom_file, fontsize_adjustment)
get_date_information(ax, cmap, dicom_file, fontsize_adjustment)
get_manufacturer_information(ax, cmap, dicom_file, fontsize_adjustment)
fig, axs = plt.subplots(2, 2, figsize=(figure_size[0] / 3 * 2, figure_size[1] / 4 * 5))
DICOM_2D_image(
axs[0, 0],
"gray",
dicom_file,
f"{dicom_file.BodyPartExamined.capitalize()} CT medical image displayed using "
"colormap 'gray'",
)
DICOM_2D_image(
axs[0, 1],
"bone",
dicom_file,
f"{dicom_file.BodyPartExamined.capitalize()} CT medical image displayed using "
"colormap 'bone'",
)
DICOM_2D_image(
axs[1, 0],
"gray",
random_dicom_file,
f"{random_dicom_file.BodyPartExamined.capitalize()} CT medical image displayed using "
"colormap 'gray'",
)
DICOM_2D_image(
axs[1, 1],
"bone",
random_dicom_file,
f"{random_dicom_file.BodyPartExamined.capitalize()} CT medical image displayed using "
"colormap 'bone'",
)
fig.suptitle(
"Visualization of CT Medical Image Display with Supplementary Information",
fontsize="x-large",
x=0.5,
y=0,
)
plt.tight_layout()
plt.show()
# The dataset used in this practice project is the MRI DICOM dataset of the head of
# a 52-year-old normal male mathematics professor
# The subject of the experiment is the author, who suffers from small vertical strabismus
# (hypertropia), which is a misalignment of the eyes, which can be seen in this dataset
# This dataset was provided by Lionheart, William R.B. in 2015, available online
# (https://zenodo.org/record/16956 or http://doi.org/10.5281/zenodo.16956),
# license type is Attribution-ShareAlike 4.0 International (CC BY-SA 4.0)
# https://creativecommons.org/licenses/by-sa/4.0/legalcode
path_to_head_mri = Path(
"../Datasets/An MRI Dicom Data Set of the Head of a Normal Male Human Aged 52/"
"ST000000/SE000001/"
)
# `Path.glob` selects the given relative pattern in the directory represented by the given path
# and yields all matching files (of any type)
# Since `Path.glob` returns a generator, it is converted to a list here for easy display
all_files = list(path_to_head_mri.glob("*"))
mri_data = []
for path in all_files:
# Read the individual DICOM files one by one in the list returned by the path generator
data = pydicom.read_file(path)
mri_data.append(data)
tabulation = Form_Generator()
tabulation.heading_printer("Reading all DICOM files that match the specified pattern")
statements = [
"""
path_to_head_mri = Path(
"../Datasets/An MRI Dicom Data Set of the Head of a Normal Male Human Aged 52/"
"ST000000/SE000001/"
)
all_files = list(path_to_head_mri.glob("*"))
mri_data = []
for path in all_files:
data = pydicom.read_file(path)
mri_data.append(data)
"""
]
tabulation.statement_generator(statements)
variables = ["path_to_head_mri", "all_files"]
values = [str(path_to_head_mri), str(reprlib_rules.repr(all_files))]
tabulation.variable_generator(variables, values, 1)
expressions = [
"all_files[0].parts",
"all_files[0].parent",
"all_files[0].name",
"all_files[0].suffixes",
"all_files[0].stem",
"len(all_files)",
"mri_data[0]",
"len(mri_data)",
]
results = [
str(all_files[0].parts),
str(all_files[0].parent),
str(all_files[0].name),
str(all_files[0].suffixes),
str(all_files[0].stem),
str(len(all_files)),
str(reprlib_rules.repr(mri_data[0])),
str(len(mri_data)),
]
tabulation.expression_generator(expressions, results, 1)
Reading all DICOM files that match the specified pattern +---------------------------------------------------------+ | Statement | +---------------------------------------------------------+ | path_to_head_mri = Path( | | "../Datasets/An MRI Dicom Data Set of the Head of a | | Normal Male Human Aged 52/" | | "ST000000/SE000001/" | | ) | | all_files = list(path_to_head_mri.glob("*")) | | | | mri_data = [] | | | | for path in all_files: | | data = pydicom.read_file(path) | | mri_data.append(data) | +---------------------------------------------------------+ +------------------+------------------------------------------+ | Variable | Value | +------------------+------------------------------------------+ | path_to_head_mri | ../Datasets/An MRI Dicom Data Set of the | | | Head of a Normal Male Human Aged | | | 52/ST000000/SE000001 | | all_files | [PosixPath('../Datasets/An MRI Dicom | | | Data Set of the Head of a Normal Male | | | Human Aged | | | 52/ST000000/SE000001/MR000000'), | | | PosixPath('../Datasets/An MRI Dicom | | | Data Set of the Head of a Normal Male | | | Human Aged | | | 52/ST000000/SE000001/MR000007'), | | | PosixPath('../Datasets/An MRI Dicom | | | Data Set of the Head of a Normal Male | | | Human Aged | | | 52/ST000000/SE000001/MR000009'), | | | PosixPath('../Datasets/An MRI Dicom | | | Data Set of the Head of a Normal Male | | | Human Aged | | | 52/ST000000/SE000001/MR000008'), | | | PosixPath('../Datasets/An MRI Dicom | | | Data Set of the Head of a Normal Male | | | Human Aged | | | 52/ST000000/SE000001/MR000006'), | | | PosixPath('../Datasets/An MRI Dicom | | | Data Set of the Head of a Normal Male | | | Human Aged | | | 52/ST000000/SE000001/MR000001'), ...] | +------------------+------------------------------------------+ +-----------------------+-------------------------------------+ | Expression | Result | +-----------------------+-------------------------------------+ | all_files[0].parts | ('..', 'Datasets', 'An MRI Dicom | | | Data Set of the Head of a Normal | | | Male Human Aged 52', 'ST000000', | | | 'SE000001', 'MR000000') | | all_files[0].parent | ../Datasets/An MRI Dicom Data Set | | | of the Head of a Normal Male Human | | | Aged 52/ST000000/SE000001 | | all_files[0].name | MR000000 | | all_files[0].suffixes | [] | | all_files[0].stem | MR000000 | | len(all_files) | 27 | | mri_data[0] | Dataset.file_meta | | | ------------------------------- | | | (0002, 0000) File Meta Information | | | Group Length UL: 214 | | | (0002, 0001) Fil...Group Length | | | UL: 131084 | | | (7fe0, 0010) Pixel Data | | | OW: Array of 131072 | | | elements | | len(mri_data) | 27 | +-----------------------+-------------------------------------+
# DICOM files may not be sorted according to their actual image positions, this can
# be verified by checking Slice Location (0020,1041)
# Slice Location (0020,1041) identifies the relative position of the image plane in millimeters
unordered_slices = {
index: float(slice.SliceLocation) for index, slice in enumerate(mri_data)
}
# For better viewing and processing of 3D data stored with multiple 2D DICOM files,
# these files must be sorted; otherwise the complete scanned image will become cluttered
# and useless
mri_data_ordered = sorted(mri_data, key=lambda slice: slice.SliceLocation)
ordered_slices = {
index: float(slice.SliceLocation) for index, slice in enumerate(mri_data_ordered)
}
tabulation = Form_Generator()
tabulation.heading_printer("Sorting of read DICOM files")
statements = [
"""
unordered_slices = {
index: float(slice.SliceLocation) for index, slice in enumerate(mri_data)
}
mri_data_ordered = sorted(mri_data, key=lambda slice: slice.SliceLocation)
ordered_slices = {
index: float(slice.SliceLocation) for index, slice in enumerate(mri_data_ordered)
}
"""
]
tabulation.statement_generator(statements)
variables = ["unordered_slices", "ordered_slices"]
values = [
str(unordered_slices),
str(ordered_slices),
]
tabulation.variable_generator(variables, values, 1)
expressions = [
"mri_data_ordered[0]",
"len(mri_data_ordered)",
"len(unordered_slices)",
"len(ordered_slices)",
]
results = [
str(reprlib_rules.repr(mri_data_ordered[0])),
str(len(mri_data_ordered)),
str(len(unordered_slices)),
str(len(ordered_slices)),
]
tabulation.expression_generator(expressions, results, 1)
Sorting of read DICOM files +-----------------------------------------------------------+ | Statement | +-----------------------------------------------------------+ | unordered_slices = { | | index: float(slice.SliceLocation) for index, slice in | | enumerate(mri_data) | | } | | | | mri_data_ordered = sorted(mri_data, key=lambda slice: | | slice.SliceLocation) | | | | ordered_slices = { | | index: float(slice.SliceLocation) for index, slice in | | enumerate(mri_data_ordered) | | } | +-----------------------------------------------------------+ +------------------+------------------------------------------+ | Variable | Value | +------------------+------------------------------------------+ | unordered_slices | {0: 0.0, 1: 41.9999963629367, 2: | | | 53.9999958207213, 3: 47.9999970362677, | | | 4: 35.9999959546749, 5: | | | 5.99999663091323, 6: 137.999998321624, | | | 7: 143.999998928727, 8: | | | 71.9999961590453, 9: 89.9999955528687, | | | 10: 83.9999967682912, 11: | | | 77.999996227574, 12: 149.999999502083, | | | 13: 131.999997780749, 14: | | | 23.9999946081714, 15: 17.9999979772582, | | | 16: 11.9999973042441, 17: | | | 29.9999952815023, 18: 107.999995419197, | | | 19: 119.999996566542, 20: | | | 95.9999960937442, 21: 65.9999961939969, | | | 22: 59.9999962290673, 23: | | | 101.999994745866, 24: 125.999997173645, | | | 25: 155.999992554172, 26: | | | 113.999995959439} | | ordered_slices | {0: 0.0, 1: 5.99999663091323, 2: | | | 11.9999973042441, 3: 17.9999979772582, | | | 4: 23.9999946081714, 5: | | | 29.9999952815023, 6: 35.9999959546749, | | | 7: 41.9999963629367, 8: | | | 47.9999970362677, 9: 53.9999958207213, | | | 10: 59.9999962290673, 11: | | | 65.9999961939969, 12: 71.9999961590453, | | | 13: 77.999996227574, 14: | | | 83.9999967682912, 15: 89.9999955528687, | | | 16: 95.9999960937442, 17: | | | 101.999994745866, 18: 107.999995419197, | | | 19: 113.999995959439, 20: | | | 119.999996566542, 21: 125.999997173645, | | | 22: 131.999997780749, 23: | | | 137.999998321624, 24: 143.999998928727, | | | 25: 149.999999502083, 26: | | | 155.999992554172} | +------------------+------------------------------------------+ +-----------------------+-------------------------------------+ | Expression | Result | +-----------------------+-------------------------------------+ | mri_data_ordered[0] | Dataset.file_meta | | | ------------------------------- | | | (0002, 0000) File Meta Information | | | Group Length UL: 214 | | | (0002, 0001) Fil...Group Length | | | UL: 131084 | | | (7fe0, 0010) Pixel Data | | | OW: Array of 131072 | | | elements | | len(mri_data_ordered) | 27 | | len(unordered_slices) | 27 | | len(ordered_slices) | 27 | +-----------------------+-------------------------------------+
# Fill the 3D array in a slice-per-slice manner
full_volume = [slice.pixel_array for slice in mri_data_ordered]
full_volume = np.array(full_volume)
tabulation = Form_Generator()
tabulation.heading_printer(
"One-time extraction of the overall 3D pixel array for all slices"
)
statements = [
"""
full_volume = [slice.pixel_array for slice in mri_data_ordered]
full_volume = np.array(full_volume)
"""
]
tabulation.statement_generator(statements)
variables = ["full_volume"]
values = [str(reprlib_rules.repr(full_volume))]
tabulation.variable_generator(variables, values)
expressions = ["full_volume.shape", "full_volume.dtype"]
results = [str(full_volume.shape), str(full_volume.dtype)]
tabulation.expression_generator(expressions, results)
One-time extraction of the overall 3D pixel array for all slices +-----------------------------------------------+ | Statement | +-----------------------------------------------+ | full_volume = [slice.pixel_array for slice in | | mri_data_ordered] | | full_volume = np.array(full_volume) | +-----------------------------------------------+ +-------------+------------------------------------+ | Variable | Value | +-------------+------------------------------------+ | full_volume | array([[[0, 0, 0, ..., 0, 0, 0], | | | [0, 0, 0, ..., 0, 0, 0], | | | [0, 0, 0, ..., 0, 0, 0], | | | ..., | | | [0,... ..., | | | [0, 0, 0, ..., 0, 0, 0], | | | [0, 0, 0, ..., 0, 0, 0], | | | [0, 0, 0, ..., 0, 0, 0]]], | | | dtype=uint16) | +-------------+------------------------------------+ +-------------------+----------------+ | Expression | Result | +-------------------+----------------+ | full_volume.shape | (27, 256, 256) | | full_volume.dtype | uint16 | +-------------------+----------------+
def sort_checker(list_1):
list_2 = list_1[1:]
if all(a <= b for a, b in zip(list_1, list_2)):
return "ascending"
elif all(a >= b for a, b in zip(list_1, list_2)):
return "descending"
elif all(a == b for a, b in zip(list_1, list_2)):
return "identical"
else:
return "unordered"
def title_adder(func):
@wraps(func)
def wrapper(inputs, ax, dict_slices, title=None, **kwargs):
is_sorted = sort_checker(list(dict_slices.values()))
if title is None:
title = f"DICOM image grid with {is_sorted} slice locations"
func(inputs, ax, dict_slices, **kwargs)
ax.set_title(
title,
loc="center",
pad=10,
)
return wrapper
def grid_annotator(func):
@wraps(func)
def wrapper(inputs, ax, dict_slices, fontsize_adjustment=0, **kwargs):
x_positions_1, x_positions_2, y_positions, cmap = func(inputs, ax, **kwargs)
for (y, x_1), key in zip(
itertools.product(y_positions, x_positions_1), dict_slices.keys()
):
text = f"No.{(key + 1):>03}"
ax.text(
x_1,
y,
text,
transform=ax.transData,
fontsize=7 + fontsize_adjustment,
fontfamily="monospace",
c=annotation_color(cmap),
)
for (y, x_2), value in zip(
itertools.product(y_positions, x_positions_2), dict_slices.values()
):
text = f"{value:>5.1f} mm"
ax.text(
x_2,
y,
text,
transform=ax.transData,
fontsize=7 + fontsize_adjustment,
fontfamily="monospace",
c=annotation_color(cmap),
)
return wrapper
@title_adder
@grid_annotator
def grid_DICOM_2D_image(inputs, ax, n=None, row_size=9, pad_value=1):
if n is None:
n = inputs.shape[0]
nrows = math.ceil(n / row_size)
cmap = "bone"
# The only types supported by `torch.Tensor` are float64, float32, float16, complex64,
# complex128, int64, int32, int16, int8, uint8 and bool, not uint16
inputs = torch.Tensor(inputs.astype(int)).unsqueeze(1)
x_positions_1 = [
math.ceil(
x * inputs[:row_size].shape[3]
+ 0.075 * inputs[:row_size].shape[3]
+ pad_value * x
)
for x in range(inputs[:row_size].shape[0])
]
x_positions_2 = [
math.ceil(
x * inputs[:row_size].shape[3]
+ 0.55 * inputs[:row_size].shape[3]
+ pad_value * x
)
for x in range(inputs[:row_size].shape[0])
]
y_positions = [
math.ceil(
y * inputs[:row_size].shape[2]
+ 0.125 * inputs[:row_size].shape[2]
+ pad_value * y
)
for y in range(nrows)
]
for i in range(nrows):
if len(inputs) > row_size:
row_images = inputs[:row_size]
inputs = inputs[row_size:]
else:
row_images = inputs[:]
# Unlike the images in the MNIST dataset, the DICOM file images do not take values
# in the range (0, 1), so the parameter `normalize` needs to be set to True in order to
# move their values proportionally to the range (0, 1)
# By querying the source code of `vision/torchvision/utils.py`, it is clear that any
# single-channel image is converted to a three-channel image, and all are represented
# as tensors, as shown below:
# > if tensor.dim() == 2: # single image H x W
# > tensor = tensor.unsqueeze(0)
# > if tensor.dim() == 3: # single image
# > if tensor.size(0) == 1: # if single-channel, convert to 3-channel
# > tensor = torch.cat((tensor, tensor, tensor), 0)
# > tensor = tensor.unsqueeze(0)
# >
# > if tensor.dim() == 4 and tensor.size(1) == 1: # single-channel images
# > tensor = torch.cat((tensor, tensor, tensor), 1)
grid_row = make_grid(
row_images, nrow=row_size, normalize=True, pad_value=pad_value
)
scalar_row = grid_row.numpy().transpose(1, 2, 0)[:, :, 0]
if i == 0:
if nrows != 1:
scalar_image = scalar_row
else:
scalar_image = np.full(
(
scalar_row.shape[0],
math.ceil(scalar_row.shape[1] / inputs.shape[0] * row_size),
),
np.nan,
)
scalar_image[:, : scalar_row.shape[1]] = scalar_row
elif scalar_image.shape[1] == scalar_row.shape[1]:
scalar_image = np.concatenate((scalar_image, scalar_row), axis=0)
else:
image_to_concat = np.full(
(scalar_image.shape[0] + scalar_row.shape[0], scalar_image.shape[1]),
np.nan,
)
image_to_concat[: scalar_image.shape[0], :] = scalar_image
image_to_concat[scalar_image.shape[0] :, : scalar_row.shape[1]] = scalar_row
scalar_image = image_to_concat
ax.imshow(scalar_image, cmap=cmap)
ax.set(xticks=[], yticks=[], frame_on=False)
ax.grid(False)
return x_positions_1, x_positions_2, y_positions, cmap
unordered_full_volume = [slice.pixel_array for slice in mri_data]
unordered_full_volume = np.array(unordered_full_volume)
descending_slices = {
index: float(slice.SliceLocation)
for index, slice in enumerate(mri_data_ordered[::-1])
}
fig, axs = plt.subplots(3, 1, figsize=(figure_size[0], figure_size[1] / 7 * 9))
grid_DICOM_2D_image(unordered_full_volume, axs[0], unordered_slices)
grid_DICOM_2D_image(full_volume, axs[1], ordered_slices)
grid_DICOM_2D_image(full_volume[::-1], axs[2], descending_slices)
fig.suptitle(
"Visual Comparison of DICOM Image Grids with Unordered and Ordered Slice Locations",
fontsize="x-large",
x=0.5,
y=0,
)
plt.tight_layout()
plt.show()
fig, axs = plt.subplots(9, 3, figsize=(
figure_size[0] / 3 * 2, figure_size[1] / 2 * 7))
for slice_counter, (i, j) in enumerate(itertools.product(range(9), range(3))):
DICOM_2D_image(
axs[i][j],
"bone",
mri_data_ordered[slice_counter],
f"DICOM image with {mri_data_ordered[slice_counter].SliceLocation:^.1f} mm "
"slice location",
fontsize_adjustment=-3,
patien_age="52",
)
fig.suptitle(
"Visualization of MRI Medical Images with Information Displayed in Order of Slice Location",
fontsize="x-large",
x=0.5,
y=0,
)
plt.tight_layout()
plt.show()
# The path object must be converted to a string in order for the SimpleITK package to process
# the path information
path_string = str(path_to_head_mri)
# For some image formats such as DICOM, the images also contain associated metadata,
# which is not loaded by the reader by default for time-saving reasons
# `ImageSeriesReader` class provides the ability to read a series of image files into a
# SimpleITK image, and once a series of images has been read, the metadata can be accessed
# directly from the reader
# `GetGDCMSeriesIDs` provides the ability to get all seriesIDs from a DICOM dataset
series_ids = sitk.ImageSeriesReader.GetGDCMSeriesIDs(path_string)
# `GetGDCMSeriesFileNames` generates a series of filenames from the directory containing the
# DICOM dataset and series ID
# Note that the resulting series of filenames (which are actually file paths) are automatically
# sorted
series_file_names = sitk.ImageSeriesReader.GetGDCMSeriesFileNames(
path_string, series_ids[0]
)
series_reader = sitk.ImageSeriesReader()
series_reader.SetFileNames(series_file_names)
# Finally, the reader is executed by calling `Execute` to get the required data
image_data = series_reader.Execute()
tabulation = Form_Generator()
tabulation.heading_printer(
"Automatic recognition and recall of series filenames for DICOM images"
)
statements = [
"""
path_string = str(path_to_head_mri)
series_ids = sitk.ImageSeriesReader.GetGDCMSeriesIDs(path_string)
series_file_names = sitk.ImageSeriesReader.GetGDCMSeriesFileNames(
path_string, series_ids[0]
)
series_reader = sitk.ImageSeriesReader()
series_reader.SetFileNames(series_file_names)
image_data = series_reader.Execute()
"""
]
tabulation.statement_generator(statements)
variables = ["series_ids", "series_file_names", "series_reader", "image_data"]
values = [
str(reprlib_rules.repr(series_ids)),
str(reprlib_rules.repr(series_file_names)),
str(reprlib_rules.repr(series_reader)),
str(reprlib_rules.repr(image_data)),
]
tabulation.variable_generator(variables, values, 1)
expressions = [
"type(series_ids)",
"len(series_ids)",
"type(series_file_names)",
"len(series_file_names)",
"len(image_data)",
"image_data.GetSize()",
]
results = [
str(type(series_ids)),
str(len(series_ids)),
str(type(series_file_names)),
str(len(series_file_names)),
str(len(image_data)),
str(image_data.GetSize()),
]
tabulation.expression_generator(expressions, results)
Automatic recognition and recall of series filenames for DICOM images +----------------------------------------------------------+ | Statement | +----------------------------------------------------------+ | path_string = str(path_to_head_mri) | | series_ids = | | sitk.ImageSeriesReader.GetGDCMSeriesIDs(path_string) | | series_file_names = | | sitk.ImageSeriesReader.GetGDCMSeriesFileNames( | | path_string, series_ids[0] | | ) | | | | series_reader = sitk.ImageSeriesReader() | | series_reader.SetFileNames(series_file_names) | | | | image_data = series_reader.Execute() | +----------------------------------------------------------+ +-------------------+-----------------------------------------+ | Variable | Value | +-------------------+-----------------------------------------+ | series_ids | ('1.3.46.67058...1413262801702',) | | series_file_names | ('../Datasets/...0001/MR000000', | | | '../Datasets/...0001/MR000001', | | | '../Datasets/...0001/MR000002', | | | '../Datasets/...0001/MR000003', | | | '../Datasets/...0001/MR000004', | | | '../Datasets/...0001/MR000005', ...) | | series_reader | <SimpleITK.SimpleITK.ImageSeriesReader; | | | proxy of <Swig Object of type | | | 'itk::simple::ImageSeriesReader *' at | | | 0x29c147690> > | | image_data | <SimpleITK.SimpleITK.Image; proxy of | | | <Swig Object of type 'std::vector< | | | itk::simple::Image >::value_type *' at | | | 0x29c147a50> > | +-------------------+-----------------------------------------+ +-------------------------+-----------------+ | Expression | Result | +-------------------------+-----------------+ | type(series_ids) | ⟨class 'tuple'⟩ | | len(series_ids) | 1 | | type(series_file_names) | ⟨class 'tuple'⟩ | | len(series_file_names) | 27 | | len(image_data) | 1769472 | | image_data.GetSize() | (256, 256, 27) | +-------------------------+-----------------+
# As shown above, the shape of the original SimpleITK image is (256, 256, 27) instead of the
# common (27, 256, 256) shape, this is due to the different order of the image dimensions,
# so it's necessary to perform the last step to convert the SimpleITK image to a NumPy array
# `GetArrayFromImage` returns a copy of the image data, which can be freely modified
# as it has no effect on the original SimpleITK image
head_mri = sitk.GetArrayFromImage(image_data)
tabulation = Form_Generator()
tabulation.heading_printer("Getting the NumPy array from the SimpleITK image")
statements = ["head_mri = sitk.GetArrayFromImage(image_data)"]
tabulation.statement_generator(statements)
variables = ["head_mri"]
values = [str(reprlib_rules.repr(head_mri))]
tabulation.variable_generator(variables, values)
expressions = [
"type(head_mri)",
"len(head_mri)",
"head_mri.shape",
]
results = [
str(type(head_mri)),
str(len(head_mri)),
str(head_mri.shape),
]
tabulation.expression_generator(expressions, results)
Getting the NumPy array from the SimpleITK image +-----------------------------------------------+ | Statement | +-----------------------------------------------+ | head_mri = sitk.GetArrayFromImage(image_data) | +-----------------------------------------------+ +----------+--------------------------------------------------+ | Variable | Value | +----------+--------------------------------------------------+ | head_mri | array([[[0, 0, 0, ..., 0, 0, 0], | | | [0, 0, 0, ..., 0, 0, 0], | | | [0, 0, 0, ..., 0, 0, 0], | | | ..., | | | [0,... ..., | | | [0, 0, 0, ..., 0, 0, 0], | | | [0, 0, 0, ..., 0, 0, 0], | | | [0, 0, 0, ..., 0, 0, 0]]], dtype=uint16) | +----------+--------------------------------------------------+ +----------------+-------------------------+ | Expression | Result | +----------------+-------------------------+ | type(head_mri) | ⟨class 'numpy.ndarray'⟩ | | len(head_mri) | 27 | | head_mri.shape | (27, 256, 256) | +----------------+-------------------------+
head_mri_directory_path = Path(
"../Datasets/An MRI Dicom Data Set of the Head of a Normal Male Human Aged 52/"
"ST000000/"
)
all_paths = sorted(list(head_mri_directory_path.glob("SE*")))
path = all_paths[-1]
new_series_ids = sitk.ImageSeriesReader.GetGDCMSeriesIDs(str(path))
new_series_file_names = sitk.ImageSeriesReader.GetGDCMSeriesFileNames(
str(path), new_series_ids[0]
)
new_series_reader = sitk.ImageSeriesReader()
new_series_reader.SetFileNames(new_series_file_names)
new_image_data = new_series_reader.Execute()
new_head_mri = sitk.GetArrayFromImage(new_image_data)
# It is easy to see that although the ordered image data can be easily obtained by SimpleITK,
# if other related metadata is needed in addition to the image data, it is still necessary
# to obtain it in the conventional manner
new_all_files = list(path.glob("*"))
new_mri_data = sorted(
[pydicom.read_file(p) for p in new_all_files],
key=lambda slice: slice.SliceLocation,
)
new_ordered_slices = {
index: float(slice.SliceLocation) for index, slice in enumerate(new_mri_data)
}
fig = plt.figure(
figsize=(figure_size[0], figure_size[1]), constrained_layout=True)
gs = gridspec.GridSpec(
nrows=2, ncols=1, figure=fig, wspace=0.08, hspace=0.04, height_ratios=[2, 3]
)
ax_1 = fig.add_subplot(gs[0, 0])
grid_DICOM_2D_image(
head_mri, ax_1, ordered_slices, title="DICOM image grid for MRI axial plane"
)
ax_2 = fig.add_subplot(gs[1, 0])
grid_DICOM_2D_image(
new_head_mri,
ax_2,
new_ordered_slices,
title="DICOM image grid for MRI coronal plane",
row_size=8,
fontsize_adjustment=1,
)
fig.suptitle(
"Visualization Comparison of Different Plane DICOM Image Grids",
fontsize="x-large",
x=0.5,
y=-0.02,
)
plt.show()
# Create a new function based on the one used earlier, just to get the independent
# plane information
def get_plane(dicom_file):
IOP = dicom_file.ImageOrientationPatient
row_xyz = IOP[:3]
column_xyz = IOP[3:]
plane = [abs(round(x)) for x in np.cross(row_xyz, column_xyz)]
if plane[0] == 1 and plane[1] == 0 and plane[2] == 0:
return "sagittal"
elif plane[0] == 0 and plane[1] == 1 and plane[2] == 0:
return "coronal"
elif plane[0] == 0 and plane[1] == 0 and plane[2] == 1:
return "axial"
else:
"Unknown"
path = all_paths[0]
all_files_0 = list(path.glob("*"))
mri_data_0 = sorted(
[pydicom.read_file(p) for p in all_files_0],
key=lambda slice: slice.SliceLocation,
)
fig = plt.figure(
figsize=(figure_size[0], figure_size[1]), constrained_layout=True)
gs = gridspec.GridSpec(
nrows=2, ncols=6, figure=fig, wspace=0.08, hspace=0.04, height_ratios=[1, 1]
)
ax_1 = plt.subplot(gs[0, :2])
ax_2 = plt.subplot(gs[0, 2:4])
ax_3 = plt.subplot(gs[0, 4:])
ax_4 = plt.subplot(gs[1, 1:3])
ax_5 = plt.subplot(gs[1, 3:5])
axs = [ax_1, ax_2, ax_3, ax_4, ax_5]
for slice_counter, ax in enumerate(axs):
DICOM_2D_image(
ax,
"bone",
mri_data_0[slice_counter],
f"MRI {get_plane(mri_data_0[slice_counter])} plane image at "
f"{mri_data_0[slice_counter].SliceLocation:^.1f} mm slice location",
fontsize_adjustment=-1,
patien_age="52",
)
fig.suptitle(
"Visualization Comparison of DICOM Images in Different Planes",
fontsize="x-large",
x=0.5,
y=-0.02,
)
plt.show()
path_to_dicom = (
"../Datasets/An MRI Dicom Data Set of the Head of a Normal Male Human Aged 52/"
"ST000000/SE000001/"
)
output_folder = (
"../Datasets/An MRI Dicom Data Set of the Head of a Normal Male Human Aged 52/"
"ST000000/"
)
# The `dicom2nifti.convert_dir` module provides the `convert_directory` function, which is
# automatically imported through the package initialization, to sort all DICOM files one-by-one
# by series, and then convert them to NIfTI format
# By default, the `convert_directory` function is set to True for the `compression` parameter,
# which is used to enable or disable Gzip compression, and for the `reorient` parameter,
# which is used to reorient the DICOM images according to the LAS direction
dicom2nifti.convert_directory(path_to_dicom, output_folder)
# NIfTI stands for Neuroimaging Informatics Technology Initiative and is co-sponsored by
# the National Institute of Mental Health and the National Institute of Neurological Disorders
# and Stroke
# NIfTI defines a neuroimaging data file format designed to meet the needs of the fMRI research
# community, and although some software continues to use other file formats, NIfTI has become
# the most widely used storage standard for fMRI and other MRI research data files
# NIfTI images can be compressed using a standard open-source algorithm called Gzip, which
# greatly reduces file size and therefore the amount of storage required for imaging data
# NIfTI defines a single image file ending in the `.nii` extension, while the compressed
# version appends the `.gz` extension, i.e. `.nii.gz`
nifti = nib.load(f"{output_folder}201_t2w_tse.nii.gz")
tabulation = Form_Generator()
tabulation.heading_printer(
"Getting and extracting the NIfTI file through conversion")
statements = [
"""
path_to_dicom = (
"../Datasets/An MRI Dicom Data Set of the Head of a Normal Male Human Aged 52/"
"ST000000/SE000001/"
)
output_folder = (
"../Datasets/An MRI Dicom Data Set of the Head of a Normal Male Human Aged 52/"
"ST000000/"
)
dicom2nifti.convert_directory(path_to_dicom, output_folder)
nifti = nib.load(f"{output_folder}201_t2w_tse.nii.gz")
"""
]
tabulation.statement_generator(statements)
variables = ["nifti"]
values = [str(nifti)]
tabulation.variable_generator(variables, values, 3)
expressions = [
"len(list(nifti.header))",
"nifti.header['qoffset_x']",
"nifti.header.get_data_shape()",
"nifti.shape",
"nifti.affine",
"nifti.affine.shape",
]
results = [
str(len(list(nifti.header))),
str(nifti.header["qoffset_x"]),
str(nifti.header.get_data_shape()),
str(nifti.shape),
str(nifti.affine),
str(nifti.affine.shape),
]
tabulation.expression_generator(expressions, results, 3)
Getting and extracting the NIfTI file through conversion +-------------------------------------------------------------+ | Statement | +-------------------------------------------------------------+ | path_to_dicom = ( | | "../Datasets/An MRI Dicom Data Set of the Head of a | | Normal Male Human Aged 52/" | | "ST000000/SE000001/" | | ) | | output_folder = ( | | "../Datasets/An MRI Dicom Data Set of the Head of a | | Normal Male Human Aged 52/" | | "ST000000/" | | ) | | dicom2nifti.convert_directory(path_to_dicom, output_folder) | | | | nifti = nib.load(f"{output_folder}201_t2w_tse.nii.gz") | +-------------------------------------------------------------+ +----------+--------------------------------------------------+ | Variable | Value | +----------+--------------------------------------------------+ | nifti | ⟨class 'nibabel.nifti1.Nifti1Image'⟩ | | | data shape (256, 256, 27) | | | affine: | | | [[-9.36898887e-01 3.20514254e-02 | | | 6.37919828e-02 1.15272324e+02] | | | [ 3.03588901e-02 9.27917123e-01 | | | -8.33337128e-01 -9.72956390e+01] | | | [ 1.43172191e-02 1.29802674e-01 | | | 5.94150448e+00 -8.23735046e+01] | | | [ 0.00000000e+00 0.00000000e+00 | | | 0.00000000e+00 1.00000000e+00]] | | | metadata: | | | ⟨class 'nibabel.nifti1.Nifti1Header'⟩ object, | | | endian='<' | | | sizeof_hdr : 348 | | | data_type : b'' | | | db_name : b'' | | | extents : 0 | | | session_error : 0 | | | regular : b'' | | | dim_info : 0 | | | dim : [ 3 256 256 27 1 1 1 | | | 1] | | | intent_p1 : 0.0 | | | intent_p2 : 0.0 | | | intent_p3 : 0.0 | | | intent_code : none | | | datatype : uint16 | | | bitpix : 16 | | | slice_start : 0 | | | pixdim : [-1. 0.9375 0.9375 | | | 5.9999995 1. 1. | | | 1. 1. ] | | | vox_offset : 0.0 | | | scl_slope : nan | | | scl_inter : nan | | | slice_end : 0 | | | slice_code : unknown | | | xyzt_units : 2 | | | cal_max : 0.0 | | | cal_min : 0.0 | | | slice_duration : 0.0 | | | toffset : 0.0 | | | glmax : 0 | | | glmin : 0 | | | descrip : b'' | | | aux_file : b'' | | | qform_code : unknown | | | sform_code : aligned | | | quatern_b : -0.016685799 | | | quatern_c : -0.9974202 | | | quatern_d : -0.069515765 | | | qoffset_x : 115.27232 | | | qoffset_y : -97.29564 | | | qoffset_z : -82.373505 | | | srow_x : [-9.3689889e-01 3.2051425e-02 | | | 6.3791983e-02 1.1527232e+02] | | | srow_y : [ 3.035889e-02 9.279171e-01 | | | -8.333371e-01 -9.729564e+01] | | | srow_z : [ 1.4317219e-02 1.2980267e-01 | | | 5.9415045e+00 -8.2373505e+01] | | | intent_name : b'' | | | magic : b'n+1' | +----------+--------------------------------------------------+ +-------------------------------+---------------------+ | Expression | Result | +-------------------------------+---------------------+ | len(list(nifti.header)) | 43 | | nifti.header['qoffset_x'] | 115.27232 | | nifti.header.get_data_shape() | (256, 256, 27) | | nifti.shape | (256, 256, 27) | | nifti.affine | [[-9.36898887e-01 | | | 3.20514254e-02 | | | 6.37919828e-02 | | | 1.15272324e+02] | | | [ 3.03588901e-02 | | | 9.27917123e-01 | | | -8.33337128e-01 | | | -9.72956390e+01] | | | [ 1.43172191e-02 | | | 1.29802674e-01 | | | 5.94150448e+00 | | | -8.23735046e+01] | | | [ 0.00000000e+00 | | | 0.00000000e+00 | | | 0.00000000e+00 | | | 1.00000000e+00]] | | nifti.affine.shape | (4, 4) | +-------------------------------+---------------------+
# Similar to DICOM files, image pixel data can be extracted from the NIfTI file using the
# `get_fdata` function
image_array = nifti.get_fdata()
tabulation = Form_Generator()
tabulation.heading_printer("Getting pixel data from the NIfTI file")
statements = ["image_array = nifti.get_fdata()"]
tabulation.statement_generator(statements)
variables = ["image_array"]
values = [str(reprlib_rules.repr(image_array))]
tabulation.variable_generator(variables, values)
expressions = [
"image_array.dtype",
"image_array.shape",
]
results = [
str(image_array.dtype),
str(image_array.shape),
]
tabulation.expression_generator(expressions, results)
Getting pixel data from the NIfTI file +---------------------------------+ | Statement | +---------------------------------+ | image_array = nifti.get_fdata() | +---------------------------------+ +-------------+------------------------------------------+ | Variable | Value | +-------------+------------------------------------------+ | image_array | array([[[0., 0., 0., ..., 0., 0., 0.], | | | [0., 0., 0., ..., 0., 0., 0.], | | | [0., 0., 0., ..., 0., 0., 0.], | | | ... ..., | | | [0., 0., 0., ..., 0., 0., 0.], | | | [0., 0., 0., ..., 0., 0., 0.], | | | [0., 0., 0., ..., 0., 0., 0.]]]) | +-------------+------------------------------------------+ +-------------------+----------------+ | Expression | Result | +-------------------+----------------+ | image_array.dtype | float64 | | image_array.shape | (256, 256, 27) | +-------------------+----------------+
def image_display_NIfTI(image, ax, title, rotation=0, cmap="bone"):
# `scipy.ndimage.rotate` uses spline interpolation of the requested order to return
# an array rotated in the plane defined by the two axes given by the `axes` parameter
# By default, the `reshape` parameter is set to True, indicating that the output shape
# is adjusted so that the input array is fully contained in the output
rotated_image = ndimage.rotate(image, rotation, reshape=True)
ax.imshow(rotated_image, cmap)
ax.grid(False)
ax.set_title(title, loc="center", pad=10)
x_ticks = list(range(0, image.shape[1], 50))
y_ticks = list(range(0, image.shape[0], 50))
ax.set(xticks=x_ticks, xticklabels=x_ticks,
yticks=y_ticks, yticklabels=y_ticks)
ax.set_xlim(left=0)
ax.set_ylim(top=0)
ax.minorticks_on()
return ax
slices_list = np.linspace(
0, image_array.shape[-1], num=4, endpoint=False, dtype=int)
plt.rcParams["figure.figsize"] = (
figure_size[0] / 3 * 2, figure_size[1] / 4 * 5)
fig, axs = plt.subplots(nrows=2, ncols=2)
for i in range(4):
if i == 0:
title = "Unrotated NIfTI single image"
else:
title = f"NIfTI single image rotated {i * 90} degrees"
image_display_NIfTI(
image_array[:, :, slices_list[i]],
axs[i // 2, i % 2],
title,
rotation=i * 90,
cmap="bone",
)
fig.suptitle(
"Visual Comparison of NIfTI Single Image at Different Rotation Angles",
fontsize="x-large",
x=0.5,
y=0,
)
plt.tight_layout()
plt.show()
def NIfTI_title_adder(func):
@wraps(func)
def wrapper(*args, rotation=0, title=None, **kwargs):
if title is None:
if rotation == 0:
title = "Unrotated NIfTI image grid"
else:
title = f"NIfTI image grid with rotated {rotation} degrees"
ax = func(*args, rotation, **kwargs)
ax.set_title(
title,
loc="center",
pad=10,
)
return wrapper
@NIfTI_title_adder
def grid_NIfTI_2D_image(
inputs,
ax,
rotation=0,
slice_axis=0,
slice_interval=1,
n=None,
row_size=9,
pad_value=1,
):
# As in the case of the SimpleITK package, NIfTI data converted from DICOM data changes
# the dimensional order of the 3D volume after conversion, resulting in a shape of
# (256, 256, 27) instead of the common (27, 256, 256) shape
# In fact, the 3D image in NIfTI data (3D or 4D) can be regarded as a stack of 2D images
# along three axes respectively, so it seems that the order of the three dimensions, i.e.,
# the three axes, is no longer so important here, but the important thing is the way to
# extract the stacked 2D images along any axis
if slice_axis == 2:
inputs = inputs.transpose(2, 0, 1)
elif slice_axis == 1:
inputs = inputs.transpose(1, 0, 2)
if n is None:
n = math.ceil(inputs.shape[0] / slice_interval)
else:
n = math.ceil(n / slice_interval)
nrows = math.ceil(n / row_size)
cmap = "bone"
rotated_inputs = [
ndimage.rotate(inputs[i, :, :], rotation, reshape=True)
for i in range(len(inputs))
]
rotated_inputs = np.array(rotated_inputs)
inputs = torch.Tensor(rotated_inputs.astype(int)).unsqueeze(1)
inputs = inputs[::slice_interval, :, :, :]
for i in range(nrows):
if len(inputs) > row_size:
row_images = inputs[:row_size]
inputs = inputs[row_size:]
else:
row_images = inputs[:]
grid_row = make_grid(
row_images, nrow=row_size, normalize=True, pad_value=pad_value
)
scalar_row = grid_row.numpy().transpose(1, 2, 0)[:, :, 0]
if i == 0:
if nrows != 1:
scalar_image = scalar_row
else:
scalar_image = np.full(
(
scalar_row.shape[0],
math.ceil(
scalar_row.shape[1] / inputs.shape[0] * row_size),
),
np.nan,
)
scalar_image[:, : scalar_row.shape[1]] = scalar_row
elif scalar_image.shape[1] == scalar_row.shape[1]:
scalar_image = np.concatenate((scalar_image, scalar_row), axis=0)
else:
image_to_concat = np.full(
(scalar_image.shape[0] +
scalar_row.shape[0], scalar_image.shape[1]),
np.nan,
)
image_to_concat[: scalar_image.shape[0], :] = scalar_image
image_to_concat[scalar_image.shape[0]:,
: scalar_row.shape[1]] = scalar_row
scalar_image = image_to_concat
ax.imshow(scalar_image, cmap=cmap)
ax.set(xticks=[], yticks=[], frame_on=False)
ax.grid(False)
return ax
fig, axs = plt.subplots(4, 1, figsize=(
figure_size[0], figure_size[1] / 7 * 12))
for i in range(4):
grid_NIfTI_2D_image(image_array, axs[i], rotation=i * 90, slice_axis=2)
fig.suptitle(
"Visual Comparison of NIfTI Image Grid at Different Rotation Angles",
fontsize="x-large",
x=0.5,
y=0,
)
plt.tight_layout()
plt.show()
# In this preprocessing step, a simple threshold is set to turn all image voxels with values
# not greater than 300 to 0
# The principle is to determine the True and False values of all image voxels and multiply them
# by the corresponding value of 1 or 0, respectively
image_array_processed = image_array * (image_array > 300)
tabulation = Form_Generator()
tabulation.heading_printer("A simple preprocessing step")
statements = ["image_array_processed = image_array * (image_array > 300)"]
tabulation.statement_generator(statements)
variables = ["image_array_processed"]
values = [str(reprlib_rules.repr(image_array_processed))]
tabulation.variable_generator(variables, values, 9)
expressions = [
"image_array > 300",
"image_array.dtype",
"image_array.shape",
]
results = [
str(reprlib_rules.repr(image_array > 300)),
str(image_array_processed.dtype),
str(image_array_processed.shape),
]
tabulation.expression_generator(expressions, results, 9)
A simple preprocessing step +-----------------------------------------------------------+ | Statement | +-----------------------------------------------------------+ | image_array_processed = image_array * (image_array > 300) | +-----------------------------------------------------------+ +-----------------------+-----------------------------------+ | Variable | Value | +-----------------------+-----------------------------------+ | image_array_processed | array([[[0., 0., 0., ..., 0., 0., | | | 0.], | | | [0., 0., 0., ..., 0., 0., | | | 0.], | | | [0., 0., 0., ..., 0., 0., | | | 0.], | | | ... ..., | | | [0., 0., 0., ..., 0., 0., | | | 0.], | | | [0., 0., 0., ..., 0., 0., | | | 0.], | | | [0., 0., 0., ..., 0., 0., | | | 0.]]]) | +-----------------------+-----------------------------------+ +-------------------+------------------------------------+ | Expression | Result | +-------------------+------------------------------------+ | image_array > 300 | array([[[False, False, False, ..., | | | False, False, False], | | | [False, False, False, ..., | | | False, False, False], | | | [... False], | | | [False, False, False, ..., | | | False, False, False], | | | [False, False, False, ..., | | | False, False, False]]]) | | image_array.dtype | float64 | | image_array.shape | (256, 256, 27) | +-------------------+------------------------------------+
plt.rcParams["figure.figsize"] = (
figure_size[0] / 3 * 2, figure_size[1] / 4 * 5)
fig, axs = plt.subplots(nrows=2, ncols=2)
for i in range(4):
if i == 0:
title = "Unrotated preprocessed NIfTI single image"
else:
title = f"Preprocessed NIfTI single image rotated {i * 90} degrees"
image_display_NIfTI(
image_array_processed[:, :, slices_list[i]],
axs[i // 2, i % 2],
title,
rotation=i * 90,
cmap="bone",
)
fig.suptitle(
"Visual Comparison of Preprocessed NIfTI Single Images with Different Rotation Angles",
fontsize="x-large",
x=0.5,
y=0,
)
plt.tight_layout()
plt.show()
# `nib.Nifti1Image` returns the initialized image, which is a combination of array-like object,
# affine matrix, header and optional metadata in `extra`, et cetera
# A NiBabel image consists of three parts: a 3D or 4D array of image data, an affine array
# that tells the position of the image array data in a reference space, and image metadata
# to describe the image
# In particular, the affine array stores a relationship between voxel coordinates in the image
# data array and coordinates in the reference space
processed_nifti = nib.Nifti1Image(image_array_processed, nifti.affine)
# `nib.save` saves NIfTI image with user-defined filename
nib.save(processed_nifti, f"{output_folder}201_t2w_tse_processed.nii.gz")
tabulation = Form_Generator()
tabulation.heading_printer("Initialization and saving of NIfTI image")
statements = [
"""
processed_nifti = nib.Nifti1Image(image_array_processed, nifti.affine)
nib.save(processed_nifti, f"{output_folder}201_t2w_tse_processed.nii.gz")
"""
]
tabulation.statement_generator(statements)
variables = ["processed_nifti"]
values = [str(processed_nifti)]
tabulation.variable_generator(variables, values, 3)
expressions = [
"len(list(processed_nifti.header))",
"processed_nifti.shape",
"processed_nifti.affine",
"processed_nifti.affine.shape",
"processed_nifti.extra",
]
results = [
str(len(list(processed_nifti.header))),
str(processed_nifti.shape),
str(processed_nifti.affine),
str(processed_nifti.affine.shape),
str(processed_nifti.extra),
]
tabulation.expression_generator(expressions, results, 3)
Initialization and saving of NIfTI image +----------------------------------------------------------+ | Statement | +----------------------------------------------------------+ | processed_nifti = nib.Nifti1Image(image_array_processed, | | nifti.affine) | | | | nib.save(processed_nifti, | | f"{output_folder}201_t2w_tse_processed.nii.gz") | +----------------------------------------------------------+ +-----------------+-------------------------------------------+ | Variable | Value | +-----------------+-------------------------------------------+ | processed_nifti | ⟨class 'nibabel.nifti1.Nifti1Image'⟩ | | | data shape (256, 256, 27) | | | affine: | | | [[-9.36898887e-01 3.20514254e-02 | | | 6.37919828e-02 1.15272324e+02] | | | [ 3.03588901e-02 9.27917123e-01 | | | -8.33337128e-01 -9.72956390e+01] | | | [ 1.43172191e-02 1.29802674e-01 | | | 5.94150448e+00 -8.23735046e+01] | | | [ 0.00000000e+00 0.00000000e+00 | | | 0.00000000e+00 1.00000000e+00]] | | | metadata: | | | ⟨class 'nibabel.nifti1.Nifti1Header'⟩ | | | object, endian='<' | | | sizeof_hdr : 348 | | | data_type : b'' | | | db_name : b'' | | | extents : 0 | | | session_error : 0 | | | regular : b'' | | | dim_info : 0 | | | dim : [ 3 256 256 27 1 | | | 1 1 1] | | | intent_p1 : 0.0 | | | intent_p2 : 0.0 | | | intent_p3 : 0.0 | | | intent_code : none | | | datatype : float64 | | | bitpix : 64 | | | slice_start : 0 | | | pixdim : [-1. | | | 0.93749994 0.9375 5.9999995 1. | | | 1. | | | 1. 1. ] | | | vox_offset : 0.0 | | | scl_slope : nan | | | scl_inter : nan | | | slice_end : 0 | | | slice_code : unknown | | | xyzt_units : 0 | | | cal_max : 0.0 | | | cal_min : 0.0 | | | slice_duration : 0.0 | | | toffset : 0.0 | | | glmax : 0 | | | glmin : 0 | | | descrip : b'' | | | aux_file : b'' | | | qform_code : unknown | | | sform_code : aligned | | | quatern_b : -0.016685799 | | | quatern_c : -0.9974202 | | | quatern_d : -0.06951577 | | | qoffset_x : 115.27232 | | | qoffset_y : -97.29564 | | | qoffset_z : -82.373505 | | | srow_x : [-9.3689889e-01 | | | 3.2051425e-02 6.3791983e-02 | | | 1.1527232e+02] | | | srow_y : [ 3.035889e-02 | | | 9.279171e-01 -8.333371e-01 | | | -9.729564e+01] | | | srow_z : [ 1.4317219e-02 | | | 1.2980267e-01 5.9415045e+00 | | | -8.2373505e+01] | | | intent_name : b'' | | | magic : b'n+1' | +-----------------+-------------------------------------------+ +-----------------------------------+---------------------+ | Expression | Result | +-----------------------------------+---------------------+ | len(list(processed_nifti.header)) | 43 | | processed_nifti.shape | (256, 256, 27) | | processed_nifti.affine | [[-9.36898887e-01 | | | 3.20514254e-02 | | | 6.37919828e-02 | | | 1.15272324e+02] | | | [ 3.03588901e-02 | | | 9.27917123e-01 | | | -8.33337128e-01 | | | -9.72956390e+01] | | | [ 1.43172191e-02 | | | 1.29802674e-01 | | | 5.94150448e+00 | | | -8.23735046e+01] | | | [ 0.00000000e+00 | | | 0.00000000e+00 | | | 0.00000000e+00 | | | 1.00000000e+00]] | | processed_nifti.affine.shape | (4, 4) | | processed_nifti.extra | {} | +-----------------------------------+---------------------+
# The dataset used in this practice project was collected as part of the Information eXtraction
# from Images (IXI) project, which collects MRI data from nearly 600 normal healthy subjects
# The MRI acquisition protocol for each subject consisted of T1, T2, and PD-weighted images,
# MRA images, and diffusion-weighted images (15 directions)
# The data in this dataset were acquired at three different hospitals in London using 1.5T and
# 3T MRI scanning systems
# This dataset was provided by IXI Dataset (https://brain-development.org/ixi-dataset/),
# license type is Attribution-ShareAlike 3.0 Unported (CC BY-SA 3.0)
# https://creativecommons.org/licenses/by-sa/3.0/legalcode
brain_mri_path = "../Datasets/IXI-T1/"
brain_mri = nib.load(f"{brain_mri_path}IXI662-Guys-1120-T1.nii.gz")
brain_mri_data = brain_mri.get_fdata()
shape = brain_mri.shape
# The affine matrix A describes the mapping from pixel coordinates to scanner or world
# coordinates and is a 4 x 4 matrix that can be represented as follows:
# A = | m_11, m_12, m_13, a |
# | m_21, m_22, m_23, b |
# | m_31, m_32, m_33, c |
# | 0, 0, 0, 1 |
# In the affine matrix A, a 3 x 3 matrix of m-series elements is responsible for rotating,
# scaling, and shearing the coordinates, denoted by M
# For a detailed description of rotation, scaling, translation, or shearing, see the following
# links: https://en.wikipedia.org/wiki/Affine_transformation#Image_transformation
# In the affine matrix A, there is also a 3 x 1 translation vector (a, b, c) that is
# responsible for the translation or offset of the affine matrix in three dimensions
# In this way, the pixel coordinates (i, j, k) are mapped to the coordinates (x, y, z) of the
# scanner or the world as follows:
# | x | = | m_11, m_12, m_13 | * | i | + | a' | = M * | i | + | a' |
# | y | | m_21, m_22, m_23 | | j | | b' | | j | | b' |
# | z | | m_31, m_32, m_33 | | k | | c' | | k | | c' |,
# where a' = a * i, b' = b * j, c' = c * k
# The addition of the extra row [0, 0, 0, 1] is a trick to put the translation part into the
# same matrix A as the matrix M, so that both parts can be applied by matrix multiplication
# To implement the above trick, an extra 1 needs to be added to the input (pixel) and output
# (scanner or world) coordinate vectors, like this:
# | x | = | m_11, m_12, m_13, a | * | i | = A * | i |
# | y | | m_21, m_22, m_23, b | | j | | j |
# | z | | m_31, m_32, m_33, c | | k | | k |
# | 1 | | 0, 0, 0, 1 | | 1 | | 1 |
# A more detailed explanation can be found at this link:
# https://nipy.org/nibabel/coordinate_systems.html#voxel-coordinates-are-in-voxel-space
affine = brain_mri.affine
tabulation = Form_Generator()
tabulation.heading_printer(
"Loading MRI data in NIfTI format and obtaining its affine matrix"
)
statements = [
"""
brain_mri_path = "../Datasets/IXI-T1/"
brain_mri = nib.load(f"{brain_mri_path}IXI662-Guys-1120-T1.nii.gz")
brain_mri_data = brain_mri.get_fdata()
shape = brain_mri.shape
affine = brain_mri.affine
"""
]
tabulation.statement_generator(statements)
variables = ["brain_mri", "brain_mri_data", "shape", "affine"]
values = [
str(brain_mri),
str(reprlib_rules.repr(brain_mri_data)),
str(shape),
str(affine),
]
tabulation.variable_generator(variables, values, 3)
expressions = [
"len(list(brain_mri.header))",
"brain_mri.header.get_zooms()",
"affine.shape",
"nib.aff2axcodes(affine)",
"affine[:3, :3]",
"affine[:3, 3:]",
"affine[3:, :]",
]
results = [
str(len(list(brain_mri.header))),
str(brain_mri.header.get_zooms()),
str(affine.shape),
str(nib.aff2axcodes(affine)),
str(affine[:3, :3]),
str(affine[:3, 3:]),
str(affine[3:, :]),
]
tabulation.expression_generator(expressions, results, 3)
Loading MRI data in NIfTI format and obtaining its affine matrix +-------------------------------------------------------------+ | Statement | +-------------------------------------------------------------+ | brain_mri_path = "../Datasets/IXI-T1/" | | brain_mri = | | nib.load(f"{brain_mri_path}IXI662-Guys-1120-T1.nii.gz") | | brain_mri_data = brain_mri.get_fdata() | | shape = brain_mri.shape | | affine = brain_mri.affine | +-------------------------------------------------------------+ +----------------+--------------------------------------------+ | Variable | Value | +----------------+--------------------------------------------+ | brain_mri | ⟨class 'nibabel.nifti1.Nifti1Image'⟩ | | | data shape (256, 256, 150) | | | affine: | | | [[ 1.89821944e-02 -2.72075552e-03 | | | 1.19975281e+00 -9.06798553e+01] | | | [-9.27821696e-01 1.32986516e-01 | | | 2.45456006e-02 1.02829445e+02] | | | [ 1.33014351e-01 9.28015888e-01 | | | 5.71511449e-11 -1.14823784e+02] | | | [ 0.00000000e+00 0.00000000e+00 | | | 0.00000000e+00 1.00000000e+00]] | | | metadata: | | | ⟨class 'nibabel.nifti1.Nifti1Header'⟩ | | | object, endian='<' | | | sizeof_hdr : 348 | | | data_type : b'' | | | db_name : b'' | | | extents : 0 | | | session_error : 0 | | | regular : b'r' | | | dim_info : 0 | | | dim : [ 3 256 256 150 1 1 | | | 1 1] | | | intent_p1 : 0.0 | | | intent_p2 : 0.0 | | | intent_p3 : 0.0 | | | intent_code : none | | | datatype : int16 | | | bitpix : 16 | | | slice_start : 0 | | | pixdim : [-1. 0.9375 | | | 0.9375 1.2000039 0. 0. | | | 0. 0. ] | | | vox_offset : 0.0 | | | scl_slope : nan | | | scl_inter : nan | | | slice_end : 0 | | | slice_code : unknown | | | xyzt_units : 10 | | | cal_max : 0.0 | | | cal_min : 0.0 | | | slice_duration : 0.0 | | | toffset : 0.0 | | | glmax : 0 | | | glmin : 0 | | | descrip : b'MR' | | | aux_file : b'' | | | qform_code : scanner | | | sform_code : scanner | | | quatern_b : 0.46861374 | | | quatern_c : -0.52952915 | | | quatern_d : -0.4576844 | | | qoffset_x : -90.679855 | | | qoffset_y : 102.829445 | | | qoffset_z : -114.823784 | | | srow_x : [ 1.8982194e-02 | | | -2.7207555e-03 1.1997528e+00 | | | -9.0679855e+01] | | | srow_y : [-9.27821696e-01 | | | 1.32986516e-01 2.45456006e-02 | | | 1.02829445e+02] | | | srow_z : [ 1.33014351e-01 | | | 9.28015888e-01 5.71511449e-11 | | | -1.14823784e+02] | | | intent_name : b'' | | | magic : b'n+1' | | brain_mri_data | array([[[0., 0., 0., ..., 0., 0., 0.], | | | [0., 0., 0., ..., 0., 0., 0.], | | | [0., 0., 0., ..., 0., 0., 0.], | | | ... ..., | | | [0., 0., 0., ..., 0., 0., 0.], | | | [0., 0., 0., ..., 0., 0., 0.], | | | [0., 0., 0., ..., 0., 0., 0.]]]) | | shape | (256, 256, 150) | | affine | [[ 1.89821944e-02 -2.72075552e-03 | | | 1.19975281e+00 -9.06798553e+01] | | | [-9.27821696e-01 1.32986516e-01 | | | 2.45456006e-02 1.02829445e+02] | | | [ 1.33014351e-01 9.28015888e-01 | | | 5.71511449e-11 -1.14823784e+02] | | | [ 0.00000000e+00 0.00000000e+00 | | | 0.00000000e+00 1.00000000e+00]] | +----------------+--------------------------------------------+ +------------------------------+-----------------------------+ | Expression | Result | +------------------------------+-----------------------------+ | len(list(brain_mri.header)) | 43 | | brain_mri.header.get_zooms() | (0.9375, 0.9375, 1.2000039) | | affine.shape | (4, 4) | | nib.aff2axcodes(affine) | ('P', 'S', 'R') | | affine[:3, :3] | [[ 1.89821944e-02 | | | -2.72075552e-03 | | | 1.19975281e+00] | | | [-9.27821696e-01 | | | 1.32986516e-01 | | | 2.45456006e-02] | | | [ 1.33014351e-01 | | | 9.28015888e-01 | | | 5.71511449e-11]] | | affine[:3, 3:] | [[ -90.67985535] | | | [ 102.82944489] | | | [-114.82378387]] | | affine[3:, :] | [[0. 0. 0. 1.]] | +------------------------------+-----------------------------+
fig = plt.figure(figsize=(figure_size[0], figure_size[1] * 2), constrained_layout=True)
gs = gridspec.GridSpec(
nrows=3, ncols=1, figure=fig, wspace=0.08, hspace=0.04, height_ratios=[50, 50, 49]
)
ax_1 = plt.subplot(gs[0, 0])
ax_2 = plt.subplot(gs[1, 0])
ax_3 = plt.subplot(gs[2, 0])
axs = [ax_1, ax_2, ax_3]
for i in range(3):
plane_dict = {
"S": "axial",
"I": "axial",
"L": "sagittal",
"R": "sagittal",
"A": "coronal",
"P": "coronal",
}
title = f"NIfTI image grid on the {plane_dict[nib.aff2axcodes(affine)[i]]} axis"
slice_interval, row_size = (8, 8) if shape[i] == 256 else (5, 6)
# The 3D image can be viewed as a stack of planes on three different axes, and each axis
# plane corresponds to a different pairwise combination of the three coordinate systems
# Since the first two axes have an overlapping coordinate system, their planes are oriented
# in the same direction, while the plane of the third axis will be rotated with respect to
# the first plane
# Although the customized function provides a way to rotate the image, it is also possible
# to get the image in the target orientation by simply swapping the axes and inverting
# the axis ordering
if i == 0:
# `np.flip` reverses the order of the elements in the array along the given axes, where
# the argument `axis` represents the axis to be flipped, and the default value is None,
# which means that all axes of the input array will be flipped
new_brain_mri = np.flip(brain_mri_data, axis=1)
new_brain_mri = np.flip(new_brain_mri, axis=2)
elif i == 1:
new_brain_mri = np.flip(brain_mri_data, 0)
new_brain_mri = np.flip(new_brain_mri, 2)
else:
# `np.swapaxes` interchanges the two axes of the array, where the arguments `axis1` and
# `axis2` represent the first axis and the second axis, respectively
brain_mri_swapped = np.swapaxes(brain_mri_data, axis1=0, axis2=1)
new_brain_mri = np.flip(brain_mri_swapped, 0)
grid_NIfTI_2D_image(
new_brain_mri,
axs[i],
slice_axis=i,
slice_interval=slice_interval,
row_size=row_size,
title=title,
)
fig.suptitle(
"Different NIfTI Image Grids on Three Axes",
fontsize="x-large",
x=0.5,
y=0,
)
plt.show()